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ABSTRACT 

The  one-speed,  steady  state,  diffusion  equation  was  solved  for  a  point 
source  and  for  a  normal  pencil  beam  in  a  homogeneous,  isotropically 
scattering  slab.  Numerical  results  obtained  using  diffusion  theory 
were  compared  to  available  transport  theory  results  for  two  slab  thick- 
nesses. 

This  comparison  demonstrated  that  the  diffusion  theory  approximation 
to  the  transport  equation  will  yield  accurate  results  except  within 
about  one  half  mean  free  path  of  a  boundary  and  except  within  about 
three  mean  free  paths  of  a  source.  The  best  agreement  between  diffusion 
theory  and  transport  theory  is  obtained  if  the  first  radial  buckling 
constant  in  the  diffusion  solution  is  chosen  equal  to  the  radial  buckling 
computed  using  transport  theory. 
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I.   INTRODUCTION 

After  being  introduced  into  a  diffusing  medium,  neutrons  suffer 
numerous  collisions  with  the  atomic  nuclei  of  the  medium.  The  collision 
process  is  quite  complicated,  since  either  scattering,  fission,  or 
absorption  can  occur.  As  a  result  of  streaming  and  repeated  collisions, 
neutrons  are  constantly  changing  their  location  in  the  medium  as  well 
as  their  velocity  (speed  and  direction).  Indeed,  some  neutrons  disappear 
(absorption)  while  others  reappear  (fission)  while  many  simply  change 
their  speed  and  direction  (scattering).  The  motion  of  any  one  neutron 
would  appear  to  be  quite  random. 

The  statistical  concept  of  a  distribution  function  is  convenient 
for  describing  this  complicated  motion.  In  essence,  we  consider  "typical" 
neutrons  and  try  to  find  the  neutron  density  throughout  a  medium  by  using 
the  principles  of  transport  theory.  Although  these  principles  are 
straightforward  and  the  exact  equation  (Boltzmann  equation)  governing 
transport  phenomena  can  easily  be  derived,  the  solution  to  this  equation 
is  usually  quite  complicated.  [1] 

Only  in  recent  years  have  methods  been  developed  to  obtain  exact 
solutions  to  the  transport  equation.  These  methods,  while  yielding  exact 
results,  are  highly  complex  and  require  considerable  computer  effort.  [2] 
For  these  reasons  many  problems  in  reactor  physics  and  shielding  design 
have  been  approached  using  the  diffusion  approximation  to  the  transport 
equation.  However,  diffusion  theory  is  expected  to  yield  accurate  results 
only  when  the  assumptions  of  Fick's  law  apply  to  the  problem.  [3] 

Quite  recently  A.  Leonard  and  G.  Garrettson  [4]  developed  techniques 
to  solve  the  neutron  transport  equation  exactly  for  the  class  of  problems 


involving  multidimensional  neutron  sources  in  a  one-dimensional  medium. 
Because  some  numerical  data  was  available  for  a  certain  subclass,  of 
these  problems,  it  was  felt  worth-while  to  develop  a  solution  to  the 
same  subclass  of  problems  using  the  diffusion  approximation  and  to 
compare  the  results  with  the  transport  solutions. 

In  this  thesis,  the  Green's  function  for  neutron  diffusion  in  a 
homogeneous,  isotropically-scattering  slab,  surrounded  by  a  non  reflecting 
medium,  was  determined  using  the  diffusion  approximation  to  the  transport 
equation.  In  Chapter  II,  an  analytical  expression  was  obtained  for  the 
neutron  density  from  a  point  source  arbitrarily  located  in  the  slab.  This 
was  integrated  to  yield  the  neutron  density  from  a  pencil  beam  normally 
incident  to  the  slab.  A  computer  program  was  written  to  evaluate  these 
expressions,  and  some  of  its  features  are  described  in  Chapter  III.  In 
Chapter  IV,  the  numerical  results  from  transport  theory  are  compared 
with  those  of  diffusion  theory. 

This  thesis  has  two  objectives.  First,  it  was  hoped  that  by  comparing 
the  numerical  results  of  diffusion  theory  with  those  of  transport  theory, 
accurate  estimates  could  be  made  of  the  region  in  which  the  diffusion 
approximation  can  be  expected  to  yield  satisfactory  results.  Second,  it 
was  hoped  a  systematic  method  could  be  found  to  choose  the  parameters 
used  in  the  diffusion  equation  to  yield  optimal  agreement  between  the 
results  of  diffusion  theory  and  transport  theory. 


II.  THE  DIFFUSION  THEORY  SOLUTION 


A.  THE  EQUATION  OF  CONTINUITY 

Under  the  assumptions  of  one-speed,  steady-state,  and  isotropic 

scattering,  the  equation  of  continuity  in  the  diffusion  approximation 

is 

Dv2$(r)  -  (za  -  vzj*(r)  +  S(r)  =  0  2.1 

—     a     t   —     — 

where  D  is  the  diffusion  coefficient;  $(r)  is  the  one  speed  neutron  flux 

as  a  function  of  position,  z,  is  the  one-group  absorption  cross  section, 

a 

in   is  the  one-group  fission  cross  section,  v  is  the  average  number  of 

neutrons  produced  per  fission,  and  S(rJ  represents  an  arbitrary  source 

distribution  function  (those  neutrons  which  have  not  yet  suffered  a 

collision).  S(_r)  will  later  be  considered  a  point  source.  Under  the 

assumptions  of  Fick's  law,  the  diffusion  coefficient  is  given  by 

D  =  l   /3z2,  where  E  is  the  one  group  scattering  cross  section  and 

1  =  1     +  Zj.   +  e  is  the  total  cross  section  [31. 
a    f    s  J 

The  physical  interpretation  of  eqn.  (2.1)  is  as  follows: 

-Dv2$(_r)d3r  represents  the  leakage  rate  out  of  a  volume  element,  d3r, 

via  streaming;  z  $(r)d3r  is  the  absorption  rate  in  d3r;vE^$(r)d3r  is 
a  t 

the  neutron  production  rate  in  d3r  from  fission;  and  S(r_)d3r  is  the 
production  rate  in  d3r  from  any  other  source. 

To  obtain  the  diffusion  approximation  to  the  transport  equation  one 
must  make  the  following  assumptions  (Fick's  Law): 

(1 

(2: 


(4 
(5 
(6 


The  medium  is  infinite 
The  medium  is  uniform 
There  are  no  neutron  sources  in  the  medium 


Scattering  is  isotropic  in  the  laboratory  coordinate  system 
The  neutron  flux  is  a  slowly  varying  function  of  position  . 
The  neutron  flux  is  not  a  function  of  time. 


The  medium  we  consider  is  a  homogeneous  slab  of  thickness  x  <  » 
with  isotropic  scattering,  and  we  solve  the  problem  under  the  assumptions 
of  one  speed  and  steady-state.  Thus  assumptions  (2),  (4),  and  (6)  are 
satisfied  but  (1),  (3),  and  (5)  are  not. 

In  order  to  obtain  the  diffusion  equation  (2.1),  it  was  necessary  to 
assume  the  diffusing  medium  to  be  infinite  in  all  directions  [3].  How- 
ever, neutrons  arriving  from  more  than  a  few  m.f.p.  (mean  free  paths) 
contribute  little  to  the  neutron  current  density  at  a  given  point,  so 
the  diffusion  equation  is  expected  to  be  reasonably  valid  except  within 
a  few  m.f.p.  of  the  boundaries. 

The  no  source  assumption,  implying  that  all  neutrons  contributing 
to  the  neutron  density  are  the  result  of  collisions,  is  certainly  not 
valid  in  the  cases  studied.  The  slabs  we  considered  contained  highly 
singular  sources:  the  point  source  and  the  normal  pencil  beam  source. 
However,  since  few  neutrons  originating  from  either  source  can  be  expect- 
ed to  survive  more  than  a  few  m.f.p.  without  suffering  a  randomizing 
collision,  diffusion  theory  should  produce  satisfactory  results  at 
distances  greater  than  a  few  m.f.p.  from  a  source. 

The  assumption  of  slowly  varying  flux  is  not  satisfied  near  the 
sources  or  boundaries.  In  the  derivation  of  Fick's  Law,  only  first 
order  terms  in  the  Taylor  expansion  were  retained  [3].  Although  second 
order  terms  will  not  contribute  to  the  current  density,  terms  containing 
third  and  higher  order  derivatives  will  make  a  contribution.  Thus,  it 
is  necessary  to  restrict  Fick's  Law  to  regions  where  the  second  derivative 
of  the  flux  does  not  change  rapidly  with  distance.  This  specifically 
excludes  regions  near  the  sources  we  consider  since  both  the  normal 
pencil  beam  source  and  the  point  source  (both  represented  using  Dirac 
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delta  functions)  are  highly  singular.  It  should  also  be  noted  that 

since  flux  tends  to  vary  rapidly  in  strongly  absorbing  media,  it  is 

also  necessary  to  restrict  Fick's  Law  to  systems  in  which  z  <<  z  . 

a     s 

Such  a  restriction  yields  an  average  number  of  secondaries  per  collision, 

c  -— 5 f , 

z      » 

close  to  one  [1]. 

From  the  preceding  discussion,  we  can  see  that  the  diffusion  theory 
approximations  should  produce  satisfactory  results  in  regions  of  the 
slab  not  too  near  boundaries  or  sources.  Hopefully,  our  numerical  results 
will  indicate  more  precisely  the  region  in  which  diffusion  theory  is 
accurate  and  the  degree  of  inaccuracy  outside  this  region. 

B.  THE  BOUNDARY  CONDITIONS  FOR  THE  STEADY  STATE  DIFFUSION  EQUATION 
As  is  normally  done  in  reactor  boundary  problems,  the  boundary 
conditions  at  the  surface  were  chosen  such  that 

1  d$  -    1  0    0 

$  Hn   "  £  CtC 

where  -r-  is  the  normal  derivative  of  the  flux  and  i   is  the  extrapolation 
length  [3].  This  is  equivalent  to  saying  that  the  flux  vanishes  at  a 
distance  i   from  the  boundary.  The  value  of  i   was  chosen  such  that  the 
solution  to  the  diffusion  equation  (2.1)  matched  as  closely  as  possible 
the  more  rigorous  solution  to  the  transport  equation  within  the  slab 
(away  from  the  boundaries  where  the  diffusion  equation  is  valid). 
For  a  planar  free  surface,  transport  theory  shows  that  i   =  0.71  X.  gives 
the  best  match,  where  x.  is  the  transport  m.f.p.  [3].  This  choice  of 
a   should  provide  good  agreement  between  transport  theory  and  diffusion 
theory  in  the  interior  of  the  slab.  However,  since  the  diffusion  equation 
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is  not  valid  near  boundaries,  a  solution  obtained  by  this  device  will 
not  give  the  correct  density  near  the  boundaries  (e.g.,  the  flux  does 
not  really  vanish  at  a  distance  i   outside  the  surface). 

An  additional  boundary  condition  is  obtained  in  a  subcritical  medium 
(c  <  1).  The  neutron  flux  decreases  toward  zero  as  the  distance  from 
the  source  increases. 

C.  '  SOLUTION  OF  THE  DIFFUSION  EQUATION 

Consider  the  homogeneous  slab  of  thickness,  x,  which  has  a  total 
cross  section,  E,  and  which  emits  c  secondaries  per  collision.  The  slab, 
infinite  in  both  transverse  directions,  is  surrounded  by  a  vacuum  or 
pure  absorber.  [See  Fig.  1.]  Under  the  assumptions  of  one-speed,  steady 
state,  and  isotropic  scattering,  the  equation  of  continuity  in  the  dif- 
fusion approximation  is  given  by  (2.1). 

Since  $(r_),  the  one-speed  neutron  flux,  is  a  function  only  of  the 
neutron  density,  n(r_),  and  the  neutron  velocity,  v  (assumed  constant), 

$(r)  =  vn(r), 
eqn.  (2.1)  can  be  written 

Dv2  n(r)  -  (za  -  vzf)n(r)  +  S^  =  o.  2.3 

-     a    t   -    v 

Where  xe(5,,x  +  £),  and  y,  z  e(-°°,°°). 

It  is  convenient  here  to  introduce  cylindrical  coordinates, 


(X»L)>  L~   (?»$)>  ancl  ?  =  vy2  +  z2>  which  are  illustrated  in  Fig.  2. 
Since  it  is  assumed  that  the  external  source 

S(x,r)  +   0,   for  xe(£;x  +  sl)  t   <j>e(0,2ir) 
r-*°° 

the  boundary  conditions  for  our  problem  are 
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FIGURE  I.  A  BARE  HOMOGENEOUS  SLAB 
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FIGURE  2.    CYLINDRICAL   COORDINATES 
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n(0,r)  =  n(t  +  2£,r)  =  0   for  f  e(0,°°),  <j>e  (0,2tt)  2.4 

and 

n(x,r)  ->  0    for  x  e(  l9t   +  £),  (j>  e(0,2ir).  2.5 

To  obtain  the  neutron  density  n(x,r)  from  any  uncoil ided  source, 
S(x,r),  it  is  convenient  to  have  the  Green's  function,  G(x,f;  x',r/ ), 
written  here  as  a  function  of  cylindrical  coordinates  [ref.  Fig.  2]. 
Then  the  formal  solution  to  (2.3)  for  any  source  S(r_)/v  is  given  by 

O        O         O 

The  Green's  function  has  the  following  properties: 

(1)  G  satisfies  the  homogeneous  differential  equation 

Dv2G  -  (zfl  -  vEf)G  =  0  2.7 

for  all 

x  eU,x  +  £),  y,  z  e(-°°,°°),  and  x  ^  x',  r_f   r ' . 

(2)  G  satisfies  the  boundary  conditions 

G(£,r;  x',r/)  =  G(t  +  £,r;  x',r')  =  0  2.8 

for  all  xc(£,t  +  a)   and  y,z,y',z'  e(-»,»);  [5]  and  G(x,'r'  x'.r')^  0. 

(3)  G  satisfies  the  proper  jump  condition  at  r  =  r ' .  In  cylindrical 

coordinates, [9] 

G(x,r;  x',r')^->  ^-  log  |r  -  r ' |   .  2.9 

r  ■>  r ' 
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Equivalently,  G  satisfies  the  diffusion  equation 

Dv2G(r',r')  -  (zg  -  vZf)G(r;r')  +  6^^-  r1)  =  0.   2.10 
with  a  point  source,  S(r)  =  <S3(_r  -  r/),  where 

r'  =  (x',y',z'),  x'efil,!  +  a),  and  y\  z'e  (-»,»). 
63  (r-r')  is  the  three-dimensional  Dirac  delta  function  given  by 

63{r_  -  _r ' )  =  6(x  -  x')  6(y  -y')  6(z  -  z')  in  cartesian  coordinates  and 
63(_r  -  r ' )  =  62(f  -  r ' )  6(x  -  x')  in  cylindrical  coordinates,  where 

62(r-  r')  --«IL^J  6(«  .  *-).  [4] 

To  solve  for  the  Green's  function  one  can  find  the  eigenfunctions  to 
the  homogeneous  form  of  equation  (2.1)  by  separation  of  variables.  Then 
G(jr;  r')  can  be  expanded  in  terms  of  these  eigenfunctions,  and  the 
expansion  coefficients  can  be  chosen  so  that  the  jump  condition  is 
satisfied  [7]. 

In  cylindrical  coordinates  equation  (2.10)  is 

[FW?F+  kl^lx"}  G<x&  *'£')-  jr(V*f>  G(x,r;x',r') 

=  "  I  «2(£  "  £')  6(*  -  x').  2.11 

Assuming  a  homogeneous  medium  and  a  point  source  located  at  x'  on  the  x 
axis,  symmetry  dictates  that  G(x>n  x'»0)  is  a  function  of  r  but  not  of 
4>.  In  general  then,  G  is  a  function  of  |£  -  jr1 1  rather  than  r  and  r_' 
separately.  Defining 

g(x,|r  -  r' |  ;  x')  e  G(x,?;  x' ,r') 

we  obtain 


15 


V7W*-W+W}    9(x.7;x-)  -Qg(x,r';x')  =  -l^p-«(x-x')  2.12 

where  we  have  defined  the  coefficient 

Q,-A_L.  2.13 

g(x,f;x')  satisfies  the  boundary  conditions  (2.8).  For  fixed  xre(0,T) 
and  r*(0,«>),  the  homogeneous  form  of  equation  (2.11),  satisfied  by  the 
eigenfunction  ij>(x,r),  is 

{fWlF  +^1  *(x,~r)  =  Q  ♦(x.r).  2.14 

To  find  the  eigenfunction  ^(x,r)  we  separate  variables 

*(x,r)  =  X(x)R(r), 

and  substitute  into  equation  (2.14).  This  leads  to  the  following 
separated  eigenfunction  problems: 

¥*     +  L2X  =  0  2.15 

dx^ 

with  boundary  conditions 

X(0)  =  X(t  +  2l)   =  0,* 

and 

£R    +idR  +  a2R  =  0  2.16 

dr2        r  dr 

with  the  boundary  condition 


*  Note:  For  simplicity,  the  solution  is  forced  to  o  at  x  =  0  and 

x  =  t  +  2s, ,  rather  than  at  x  .=  -  i   and  x  =  t  +  i.     For  numerical  work 

and  comparison  with  transport  theory,  a  change  of  variables,  x  =  x  -  l, 
is  introduced  where  x  e(0,x)  is  inside  the  slab. 

16 


R(r)  ^  0 


where 

a2  E  Q  +  L2. 


Equation  (2.15)  has  solutions  of  the  form 

Xm(x)  =  C1  sin(Lmx)  +  C2  cos  (l^x) 
where  the  boundary  conditions  require  that  C2  =  0  and  that  the  eigen- 


values are 


L  =  -*W»   m  =  1>2'3'  "" 
m   t  +  c.i 


By  changing  variables  such  that  p  =  <y\     where  am  =  J  Q  +  L2 
equation  (2.16)  can  be  reduced  to  the  modified  Bessel's  equation 

p2d!R(pi    +pdRM    _p2     R(p)  =  0>  2.17 

Since  the  solutions  to  equation  (2.17)  are  the  zeroth  order  modified 
Bessel  functions  IQ(p)  and  KQ(p),  [6]  the  general  solution  is 

R(p)  =  C3I0(p)  +C4Ko(p).  2.18 

The  boundary  condition  requires  that  C3  =  0  since  I  (p)  —>    °°  . 

Expanding  the  Green's  function  in  terms  of  the  eigenfunctions 

Xm(x),sin(^J      and  Rj?)  h  W) 

we  obtain 

g(x,|r-r'|;  x')  -  j^  \,(x')  Xj?)   I^CIE  -  r'|).    2-19 

Since 

Vx>  y?)  -  *m(x,r) 
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satisfies  the  homogeneous  equation  (2.14)  and  the  boundary  conditions 
(2.8),  the  expansion  (2.19)  for  G(x,r_;  x',r.')  certainly  satisfies  the 
conditions  (2.7)  and  (2.8)  required  of  the  Green's  function.  It 
remains  only  to  be  shown  that  (2.19)  satisfies  the  jump  condition  (2.9). 
This  is  obtained  immediately  since  [6] 

K  (air  -  r1 1)   ■*  -  Jin  |r  -  r '  I 
ov  m»-  -  "  ^  '-  -  ' 

so  that  Rj,(r)  satisfies  the  following  differential   equation     [8]: 


*%&   .   1     dRm^     _  a2     R  m  __  .  hm 


~3F~ 


r       dr 


m    m 


2.20 


Thus  equation  (2.19)  is  the  proper  form  for  the  solution. 

Substitution  (2.19)  into  equation  (2.10),  and  using  the  Fourier  sine 
representation  [8], 


co 


«<*-*'>=    £7Ar«MrSr*)«M?$&)-       2-21 


M--  1 


as  well  as  equation  (2.20),  yields 


Vx'>=    mrhi)  sin(7^?J- 


where  we  have  utilized  the  orthogonality  property 


sin 


m-rrX 


+  2i 


sin 


nux 


+  2i 


J 


dx  = 


0  for  m  f   n 


+  2i 


for  m  =  n 


Therefore,  recalling  that  G(x,_r;  x',71 )  =  g(x,  \r_  -   ?'  |  ;xl) , 
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CO 


g(x,|r-r'|;x')  =  /^-DT32ir  sin  [7wr]    sin  (^j-)  KQ(aJr-L' | ) 

At-  I 

2.22 

is  the  required  Green's  function  for  the  point  source  problem. 

It  would,  of  course,  be  possible  to  develop  the  Green's  function 
for  the  normal  beam  problems  in  a  like  manner.  However,  this  is  not 
necessary  since  equations  (2.22)  and  (2.6)  represent  the  formal  solution 
to  (2.3)  for  any  source.  In  particular, 

p(x,?;x')  =   I    g(x,|r-r'|;  x')  62(r)e  ~zx'  dx'  f'dr'dcf,' 

o   o    I  2-23 

is  the  neutron  density  in  a  slab  from  a  pencil  beam  normally  incident  to 
the  slab  at  x  =  y  =  z  =  0,  and  this  function  satisfies  equation  (2.3) 
with  _JL  =  62(ge  -Zx*_  Thus  p  (x,  |r-f '  I  »x* )  is  the  Green's  function 
for  the  class  of  problems  involving  neutron  beams  normally  incident  to 
a  slab.  Substitution  of  (2.22)  into  (2.23)  and  subsequent  integration 
yields 

co 
jyc;.yi\.  V^m  11-(:l)mexp[-E(T+2£)]  1  .  /m-rrx  \  y  ,    yx   ?  ?. 
P(x,r,x  )  -  l^n    UmtrP  +  zU+Zz)         f  Sin(^+2TJ  Vamr}-  2*24 

To  facilitate  comparison  of  our  results  with  those  of  transport 
theory  [1],  we  measure  length  in  dimensionless  units  of  mean  free  path 
(1  m.f.p.  =  1/z)  and  shift  coordinates  in  x  and  x1,  x  ■*■  x  -  A,  so  that 
x  e(0,r)  lies  inside  the  slab.  Then  the  solutions  (2.22)  and  (2.24) 
take  the  form 

g(x,r;x-)  -  C  £  sin  (m-^±i>)  «*{*$&■)   W)     2.25 
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and 


CO 


"•(#)  Vv>.     *■* 


where  all  lengths  are  expressed  in  units  of  mean  free  path.  For  example 


tf-vf^SS-T  (rf)- 


2.27 


where  zr  =  f/(l/z)  is  the  radial  distance  expressed  in  mean  free  paths. 
For  equations  (2.26)-(2.27) ,  and  from  this  point  on,  all  distances  are 
understood  to  be  expressed  in  units  of  mean  free  path  (e.g.,  e(t+20  ■* 
t+2i    m.f.p.).  The  constants  C  and  C  will  be  adjusted  to  normalize  the 
diffusion  theory  results  to  transport  theory  results  at  some  chosen  point, 
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III.  NUMERICAL  TECHNIQUES 

A  Fortran  IV  program  (appendix  2)  was  written  to  compute  numerical 
results  using  equations  (2.25)  and  (2.26)  for  comparison  with  numerical 
data  obtained  by  transport  theory.  Since  both  equations  involve  an 
infinite  series,  it  was  necessary  to  truncate  the  series  at  some  point 
in  order  to  obtain  a  solution.  The  purpose  of  this  chapter  is  to 
justify  the  criteria  used  to  terminate  the  series  and  also  to  explain 
the  choice  of  the  particular  values  used  for  the  parameter  Q/£2. 

A.  TRUNCATION  OF  THE  SERIES 

Both  equation  (2.25)  and  (2.26)  are  infinite  series  solutions.  In 
order  to  obtain  numerical  output,  it  was  necessary  to  truncate  these 
series  at  some  point  in  the  summation  process.  That  is, 

CO  i\/-' 

where  we  choose  N  large  enough  so  that  Rm/Sn  is  less  than  some  chosen 
e,  where 


CO 


rn=  z.ymwKo<v>>  3-2 

is  the  remainder  term.  The  maximum  absolute  value  of  the  A  X  (x)  terms 

m  m 

in  both  equation  (2.25)  and  (2.26)  is  unity.  Thus  the  K  Bessel  function 

term  dominates  the  series  for  a  r  >>  1. 

m 

Since  _ 

co 

lyjn.      \<   2]  Ko(v}  3-3 

factoring  out    Maw^)   dnd  noting  that  [6] 

p->oo 
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with 


yields 


00  ■         /vr 


This  power  series  can  be  summed  to  yield  a  convenient  estimate  for  an 
upper  bound  of  the  remainder: 

KQ(gNr) 

N  ~  /ttF  ,  *  3.5 

1  -  e  ~{t+2i> 

In  the  computer  program  the  series  is  truncated  at  N  terms  with  N 

large  enough  so  that 

— ^ <  e,  3.0 

where  we  chose  e  sufficiently  small   to  give  the  desired  accuracy.     Then 
(3.5)  yields  the  upper  bound 

RN  e 

SN        1   -  exp[-gy 

for  the  remaining  fraction.  This  upper  bound  is  of  the  order  of  e  except 

when  r  <<  ,  so  that  the  choice  e  =  0.0005  should  give  at  least 

three  place  accuracy  except  for  small  7.  However,  (3.7)  represents  a 
very  conservative  estimate  of  the  upper  bound,  so  we  would  expect  good 
accuracy  with  this  choice  of  e  even  when  r  is  small.  In  fact,  for  r  =  0.05 
and  N  large  enough  so  that 

K0(aNr) 

■  °  q     <  0.0005, 
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RN 
numerical  results  indicate  that   J-*-  <  0.001,  while  (3.7) yields  the 

R      N 

N 
conservative  upper  bound   ^—  <  0.01. 

Therefore,  on  the  basis  of  this  discussion,  it  would  appear  that 

(3.6),  with  e  =  0.0005,  is  a  sufficient  truncation  requirement  to  insure 

three-place  accuracy. 

B.  CHOICE  OF  THE  PARAMETER  Q/z2 

In  equations (2.25)  -  (2.27)  the  parameter  Q/z2  appears,  where  Q  is 
defined  in  eqn.  (2.13).  This  parameter  is  a  function  of  the  one  group 
material  cross  section  of  the  slab.  However,  the  one-group  cross  sections 
depend  specifically  upon  the  energy-dependent  flux,  which  is  not  avail- 
able. Therefore,  one  would  like  to  choose  Q/z2  so  that  diffusion  theory 
will  yield  results  that  agree  as  closely  as  possible  with  the  results  of 
one-speed  transport  theory.  In  this  thesis  it  is  demonstrated  how  Q/z2 
might  be  chosen  to  accomplish  this. 

For  comparison,  this  parameter  was  chosen  in  two  ways.  The  first 
choice  assumes  diffusion  theory  for  a  medium  in  which  there  is  no  fission. 
For  the  second  choice,  it  is  noted  that  as  r  increases,  the  transport 
theory  solutions  and  the  diffusion  theory  solutions  decay  like  K  (s  r) 
and  K  (a  r),  respectively,  where  3  are  the  transport  theory  radial 
buckling  modes  [2].  To  force  diffusion  theory  to  behave  like  transport 
theory  for  large  r,  the  first  buckling  modes  in  diffusion  theory  and 
transport  theory  are  chosen  equal ,  a,   =  B-, . 

To  obtain  the  expression  for  Q/z2  in  the  first  case,  recall  the 

parameters 

z  -  vZf 
Q  ■  a  D  f  3.8 

and 
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Z   +  vZ, 
-S-s 3.9 


defined  in  Chapter  II.  In  the  diffusion  approximation  [3]  D  =  z s/3z2  , 
and  if  there  is  no  fission  (zf  =  D).  Q/z2  reduces  to 


Q  _  3za  .  3(1-c) 

z* "  z^  c 

s 


3.10 


since  c  =  £_/(£_  +  z)   in  this  case, 
s   a    s 

Recall  that  the  diffusion  approximation  demands  a  highly  scattering 
medium,  c  c^l.O,  so  c  =  0.9  is  chosen  for  the  numerical  work  presented 
here.  For  this  value  of  c,  eqn.  (3.10)  yields  the  value  Q/z2  =  0.333. 

For  the  second  case  the  parameter  Q/£2  is  obtained  by  comparing  the 
transport  theory  and  diffusion  theory  solutions.  For  large  transverse 
distances  from  the  source,  |r.-r'  |  >>  1 ,  a  relatively  simple  transport 
theory  solution  is  obtained  [2]  since  the  finite  (M  <  ■»)  point  spectrum 
modes  dominate  the  continuous  spectrum  modes: 


Gtr(r;r-)  -^  W^£  r Jx)  r-(x.)  +  0  [  .-IEt'1]    3.11 

m- 1 
where  r  (x),  m  =  1,2,---,M,  are  the  transport  eigenfunctions  for  the 

discrete  eigenvalues  0  <  3-,  <  3o  <  *  *  *  <  3M  £  1  [2].  Using  eqn.  (3.4) 

it  can  be  seen  that  the  first  mode  dominates  (3.11)  for  large  radial 

arguments,  |r-r' |  >>  1 : 

-e2lL-L'l 

Gtr(r;r')  =  J- K 0(eJr-r'|)  r,(x)  i\(x')  +  0  [  £- ri77- I . 

1  I  L~L\ 

Likewise,  the  first  mode  m  =  1  dominates  the  diffusion  theory  solution 
(2.25)  for  |r_-r/ 1  >>1 .  Therefore,  if  we  choose  the  first  radial  buckling 
modes  to  be  equal,  a,  =  £, ,  the  two  solutions  will  behave  similarly  for 

|r  -  r'  j  >>  1 . 
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With,  this  choice  for  a,  =  I  ^r  +  (~"^+7j 
expression  is  obtained  for  Q/£2: 


1/2 


,  the  following 


Sr 


-.6f- 


:+2l 


3.12 


Numerically  this  corresponds  to  Q/z2  =  0.2757,  for  c  =  0.9  and  t  =  10 
and  to  Q/e2  =  0.1918,  for  c  =  0.9  and  x  =  3.  [2] 
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TV.  DISCUSSION  OF  RESULTS  AND  CONCLUSIONS 

The  following  results  were  obtained  using  the  computer  program  discus- 
sed in  Chapter  III.  The  cases  presented  are  those  for  which  the  exact 
solution  has  been  obtained  by  G.  Garrettson  [2]  using  transport  theory. 
The  purpose  of  this  chapter  is  to  compare  the  results  obtained  using  dif- 
fusion theory  with  those  obtained  using  transport  theory  and  to  evaluate 
the  performance  of  diffusion  theory  for  problems  of  this  class.  These 
comparisons  were  done  for  cases  of  homogeneous  slabs  of  thickness  t  =  10 
and  t  =  3  m.f.p.  (mean  free  paths)  surrounded  by  a  non-reflecting  medium, 
under  the  assumption  of  one  speed,  steady  state,  and  isotropic  scatter- 
ing. All  graphical  results  (Figs.  3  to  17)  discussed  in  this  chapter 
are  found  in  Appendix  1. 

The  first  results  presented,  Figs.  3-6,  are  those  for  a  pencil  beam, 
normally  incident  at  (0,0,0)  to  a  slab  of  thickness  x  =  10  m.f.p.,  with 
a  multiplication  constant  c  =  0.9.  The  data  has  been  normalized  to  the 
transport  theory  results  at  r  =  4.4448,  x  =  5.1282  in  order  to  facilitate 
the  comparison.  Figures  7-10  depict  the  results  obtained  from  a  pencil 
beam  normally  incident  at  (0,0,0)  to  a  slab  of  thickness  t  =  3  m.f.p. 
with  c  =  0.9.  Figures  7  and  8  were  normalized  at  r  =  2.01,  x  =  1.5, 
while  Figs.  9  and  10  were  normalized  at  f  =  0.2,  x  =  1.5.  Finally,  Figs. 
11-16  represent  a  comparison  of  the  point  source  solutions  in  a  slab  of 
thickness  x  =  10,  with  c  =  0.9,  for  point  sources  located  at  depths  of 
1  and  then  5  m.f.p. 

In  all  but  Figs.  11-13  the  diffusion  solutions  are  presented  for  both 
values  of  Q/z2  discussed  in  Chapter  III.     Tn  this  chapter,  and  in  Appendix 
1,  Q/z2  is  referred  to  simply  as  Q,  since  E  =  1  in  units  of  mean  free  path. 
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A.  COMPARISON  OF  NORMAL  BEAM  SOLUTIONS. 

It  is  interesting  to  note  the  evolution  of  the  neutron  density  from 
a  pencil  beam,  N(x,r),  vs.  x  e[0,t],  as  the  radial  distance,  r,  increases. 
For  very  small  r  (ref.  Figs.  3,  9,  and  10),  the  neutron  density  tends  to 
have  the  same  shape  as  the  once-collided  neutron  density,  which  in  turn 
resembles  the  uncoil ided  density  (a  step  function  times  an  exponential). 
Transport  theory  predicts  a  wery   sharp  discontinuity  near  the  boundaries 
which  is  physically  explained  by  the  surface  leakage  [2].  Although  dif- 
fusion theory  correctly  yields  an  exponentially  decaying  density  within 
the  slab,*  it  does  not  yield  the  steep  gradient  given  by  transport  theory 
within  about  one  half  m.f.p.  of  either  surface.  The  results  of  the  two 
theories  diverge  partly  because  transport  theory  accounts  for  the 
preferential  streaming  toward  the  boundaries  while  diffusion  theory  does 
not.  In  fact,  as  previously  discussed,  the  diffusion  equation  is  simply 
not  valid  near  boundaries.  (Ref.  Chapter  II). 

For  any  x  e(0,-r),  the  diffusion  theory  solution  appears  to  be  most 
inaccurate  for  small  radial  distances,  r  <  1.  This  inaccuracy  is  apparent 
from  both  the  normal  beam  and  point  source  data  (ref.  Figs.  3,  9,  10,  11, 
and  14),  and  it  arises  partly  because  diffusion  theory  does  not  consider 
the  preferential  streaming  from  a  singular  source  in  the  region  of  the 
source.  Since  Fick's  Law  is  derived  under  the  assumption  that  there  are 
no  sources  in  the  medium  (ref.  Chapter  II)  we  should  not  expect  the 
diffusion  results  to  be  accurate  near  sources.  In  fact,  within  one  m.f.p. 
of  the  source,  the  errors  induced  by  using  diffusion  theory  may  be  as 
large  as  80  percent.  However,  beyond  about  three  m.f.p.  the  two  theories 
yield  nearly  identical  results  (if  Q  is  computed  using  eqn.  (3.12)). 


*  Recall  that  a  straight  line  on  a  semilog  plot  corresponds  to  an 

exponential . 
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B.  COMPARISON  OF  POINT  SOURCE  SOLUTIONS 

For  solutions  to  the  point  source  problem,  the  most  noticeable 
discrepancy  betv/een  the  transport  and  diffusion  data  appears  in  the 
region  within  about  one  m.f.p.  of  the  source  (ref.  Figs.  11  and  14). 
As  previously  mentioned,  diffusion  theory's  inaccuracy  near  sources 
is  partially  due  to  its  failure  to  consider  the  preferential  streaming 
in  the  region  of  the  source.  As  the  distance  from  the  source  is 
increased,  the  relative  error  decreases  unless  a  boundary  region  is 
encountered.  As  expected,  within  one  half  m.f.p.  of  the  boundary,  the 
diffusion  theory  solution  begins  to  diverge  until  a  discrepancy  of 
about  20  percent  is  noted  at  the  edges. 

C.  CHOICE  OF  THE  PARAMETER  Q 

It  should  be  noted  (ref.  Figs.  5,  7,  8,  and  16)  that  for  large 
radial  distances,  r,  the  agreement  between  transport  and  diffusion 
theory  can  be  greatly  enhanced  by  the  proper  choice  of  the  parameter 
Q  (ref.  Chapter  III).  The  diffusion  theory  yields  Q  =  3^"c^  for  a 
medium  in  which  there  is  no  fission.  However,  it  was  expected  that 


« ■ B? "  MrA 


where  e,  is  the  first  transport  eigenvalue  [2],  would  yield  better 
results.  By  examining  both  the  normal  beam  and  the  point  source  data 
for  the  two  choices  of  Q,  the  most  optimum  choice  is  apparently  given 
by  equation  (3.1 2j.  This  choice  forces  the  diffusion  solution  to  decay 
in  r  at  the  same  rate  as  the  transport  solution  (for  large  r)  because 
the  first  modes  dominate  for  f  >>1.  The  numerical  data  bears  this  out 

2 

since  the  Q  =  e?  -!   Jj  J   solution  is  nearly  identical  to  the  transport 
solution  for  t  >  3.  Unfortunately,  this  choice  of  Q  does  not  consistently 
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improve  the  results  near  the  source  (or  near  the  boundaries)  since  the 
assumptions  for  Fick's  law  are  not  valid  in  these  regions  and  since  the 
first  mode  does  not  dominate  the  solution  for  small  r..  A  plot  of  e,  vs  t 
is  given  in  G.  Garrettson's  thesis  for  c  =  0.9,  and  data  for  other  c's 
is  contained  in  an  appendix  [2]. 

A  comparison  of  results  in  the  x  =  10  and  t  =  3  slabs  (refs.  Figs. 
3-10),  indicates  that  diffusion  theory  is  more  accurate  in  thick  slabs. 
This  is  to  be  expected  since  Fick's  law  is  not  valid  near  the  boundaries, 
and  in  a  thin  slab  the  neutron  is  never  very   far  from  one  of  the  two 
boundaries. 

D.  CONCLUSIONS 

In  general,  one  can  say  that  the  neutron  density  computed  using 
diffusion  theory  can  be  expected  to  yield  reasonably  accurate  results 
except  as  indicated  in  the  shaded  areas  of  Fig.  17.  If  one  is  not 
interested  in  the  density  within  about  three  m.f.p.  from  the  source  and 
within  about  one-half  m.f.p.  from  the  edges,  diffusion  theory  can  success- 
fully be  used.  However,  if  one  is  interested  in  the  density  within  these 
regions  (e.g.,  the  shaded  areas  in  Fig.  17),  it  will  be  necessary  to 
resort  to  the  more  rigorous  and  time-consuming  techniques  of  transport 
theory  to  obtain  accurate  results. 

One  shouod  also  consider  that  the  diffusion  theory  solution  requires 
very  long  to  compute  for  small  radial  distances  from  the  source,  since 
the  smaller  )r-r' |  is,  the  more  terms  the  truncated  series  includes.  For 
example,  computer  runs  which  included  r  =  0.01  required  as  long  as  18.5 
minutes  to  complete,  which  is  comparable  to  the  computation  time  required 
for  the  most  difficult  transport  theory  solutions.  These  runs,  even 
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though  lengthy,  did  not  provide  accurate  data,  and  one  must  conclude 
that  only  the  transport  theory  solution  should  be  used  to  compute 
densities  at  wery   small  radial  distances. 

Lastly,  the  choice  of  the  parameter  Q  is  an  important  factor  in 
obtaining  agreement  between  the  results  of  diffusion  theory  and  transport 
theory  in  the  regions  where  the  diffusion  equation  is  valid.  The  choice 
(3.12)  significantly  decreases  the  discrepancy  between  the  two  results- 
at  large  radial  distances  from  the  source,  e.g.,  r  >  3  m.f.p. 

In  shielding  problems  with  beam  sources,  the  effect  of  these 
discrepancies  will  probably  be  most  noted  when  calculating  the  transmis- 
sion or  reflection  from  the  region  of  the  beam.  In  these  regions  the 
diffusion  results  are  most  inaccurate  since  the  emerging  current  from 
a  region  of  the  slab  is  a  function  of  the  neutron  density  in  that  region 
[3].  Diffusion  theory  should  give  satisfactory  results  for  both  trans- 
mission and  reflection  from  regions  not  too  near  the  beam. 

E.  SUGGESTIONS  FOR  POSSIBLE  EXTENSIONS 

For  slabs  sufficiently  thin,  no  discrete  poles  exist  for  the  transport 
equation  [2].  In  these  cases,  to  choose  a  proper  value  for  Q,  an  effec- 
tive mode  of  the  transport  equation  would  have  to  be  found  empirically. 
This  could  be  done  by  making  a  semilog  plot  of  the  transport  solution 
for  fixed  x,  as  a  function  of  r  (for  large  r) ,  and  using  the  slope  to 
determine  the  effective  decay  constant.  Mathematically,  this  would 
correspond  to  the  smallest  "effective"  pole  in  the  continuous  spectrum 
of  the  transport  operator.  These  effective  poles  should  be  investigated, 
not  only  empirically,  as  mentioned  above,  but  also  mathematically. 
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Further,  it  would  be  interesting  to  investigate  the  results  obtained 

over  a  wide  range  of  values  of  both  t  and  c.  One  would  expect  diffusion 

theory  to  be  less  accurate  for  small  t,  because  of  leakage,  and  also  for 

E   +  vEf 
small  c  =  -=— = ,  which  corresponds  to  a  highly  absorbing  slab.  The 

diffusion  approximation  is  not  expected  to  be  valid  unless  I.     <<   z  . 

a    s 

(Ref.  Chapter  II). 
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APPENDIX  1  -  Numerical  Results 

Appendix  1  contains  graphical  displays  comparing  numerical  data 
obtained  using  transport  theory  with  that  obtained  using  diffusion 
theory.  Figures  3-10  represent  semi-log  plots  of  the  neutron  density 
in  a  slab,  as  a  function  of  normal  depth,  from  a  pencil  beam  normally 
incident  to  the  slab,  for  various  radial  distances  from  the  beam  path, 
Figures  11-16  represent  semi-log  plots  of  the  neutron  density  in  a 
slab,  as  a  function  of  normal  depth,  from  a  point  source  within  the 
slab,  for  various  radial  distances  from  the  source. 

This  data  is  discussed  in  Chapter  IV. 
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APPENDIX  2  -  Fortran  IV  Program  Listing 

C       STUDY    OF    NEUTRCN    DIFFUSION 
C 

REAL*4    NUB     ,    NVAL,     NUM,     MQOB 

DIMENSION    P(1C0,1C0)     ,     PT(100,1C0) 

DIMENSION    XMESH     (100)     ,     R^SH    (100) 

RE AD (5, 101 ) EPS, X PR  I ,TAU, C, D 

WRITE(6,101)EPS,XPRI,TAU,C,D 

READ(5,111)NVAL,PTVAL,P,Z,Y 

WRITE (6, 111 )NVAL,PTVAL,R, Z,Y 
111  F0RMAT<5F10o6) 

NR  =  13 

NX  =  35 
C    THE    35    XMESH    PTS    ARE 

READ(5, lfCn  )  (RMESHC K)  ,K=1,NR) 

WRITE(6tlOCC)(QMPSH(K),K=lTN*) 

READ     (5,121)     (XM=SH(J),J=1,NX) 

-WRITF(*,121)      (XMESH<J),J=1,NXI 
121     FORMAT     (6F1C5) 
101     FORMAT( 5F1C>5 ) 
1C00     FORMAT     (5FlCo5) 

M  =  0 

TEMP=0,0 

NUM=0„C 

MODB=0,0 

PI  =  ?ol41.593 

T=TAU-M2<,*EPS) 

AR=0,0 

DCY=0o0 

ALPMR=Oa0 

FUD=0o0 

CH=030 

A  =  OoO 

XSYN=OoO 

BESL=OoO 

TO=OoO 

NUB=0,0 

CK=OoO 

C2=OoO 

Z=7+EPS 

Y=Y    +    EPS 

XP=0„0 

XP=XPRI+EPS 
C 

C    COMPUTE    THE     NORMALIZING    CONSTANTS 
5    M=M+1 

AB  =  M 

DKY=AB*( l.-( (  i-lo  )**M)*EXP(-T     )  )  )/(  ( AB*PI  1*^2  +  11**2)  ) 

C2=A8*PI/T 

A=R*SQRT(C+(C2**2)  ) 

XSYN=SIN(C2*7 t 

CALL    RESM  A,C,3K,  IFR) 

IF     (  IEWoFQoC)     GO    TO    200 

WRITE(6,201)      IER 
201    FORMATC «0« , ******    THE    BESSEL    ERROR    I S • , 13, •*****• ) 
200    CONTINUE 

□  r  c  I    =  Q  |( 

TO=T0+<DKY*XSYN*BESL) 

NUB  =  NUB  «-<2*/T)*$IN<C2*Y)*XSYN*BESL 
C  COMPUTE  SIZE  OF  BESSEL  FN 

CK=BESL/TO 

IF( ABS(CK)-Co000  5)125,5,5 
125  CONTINUE 

CCNN=NVAL/TO 

CCNP=PTVAL/NUB 

M  =  0 

DO  3  K=1,NR 

AR=RMESH(K) 

DO  4  J=1,NX 

X=XMESH( J) 

AX=X+EPS 
120  M=M+1 

AM  =  M 
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c 
c 
c 
c 

c 
c 

c 
c 
c 

c 
c 

c 
c 


CCMPUTE  DECAY  TERM,  DCY 

DCY=AM*( lo-( ( I -la ) **M)*EXPl-T) I )/ H AM+PI |**2+(T**2)) 

CCMPUTW  FRACTICN  M*PI/T=C1 
C1=AM*PI/T 


COMPUTE  ALPMR 

ALPMR=AR*SORT( (D)+(C1**2)1 

SIN  X  TER^ 

FUD=  SIN(C1*AX) 

BESSEL    FN    TERM 

CALL     RCSK( ALPMR,Ot BK, I EP ) 
IF     (  IERoEQ.>n     GO    Tn    202 
WRITE(6,201)     IER 
202    CCNTINUE 
MODB=BK 

TEMP=TEMP    f (DCY*M0DB*FUD1 

CCMPLTE  VALUE  FROM  POINT  SOURCE 

NUM=NUM+<(?0/TH:SIN(Cl*XP)  *FUD*MODB) 
CHECK  VALUE  OF  BESSEL  FN 

CH=MODR/TEMP 

IF  (ABS(CH)-O. 0005)140, 120, 120 

140  CCNTINUE 

B( J,K)=CONN*TEMP 
PT(  J,K)  =CONP*NUM 
M  =  0 

NUM=OoO 
TEMP=OoOO 
4  CONT  INUE 
3  CONTINUE 


500  FORMAT( 
WRITF(6 


501  FORMAT( //«0«  ,  iTHF* ♦  13,  •  X  MESH  POINTS  ARE'/) 


502  FORMAT( 
WRITE(6 


WRITF  (6,502) (XMESHC J) ,J=1,NX) 


503  FOPMAT( 
WRITE(6 

504  FORMAL 
WRITE(6 

505  FORMATt 


DO  105  K=1,NR 


WRITF(6 

WRITF( 6 
508  FORMAT 
105  CONTINU 

WRITE( 6 
506  FORMAT( 

WRITE  ( 6,506  ) 

DO  106  K=1,NR 

WRITE(6 


WRITE  (6,507) 


507  FORMAT 
WRITE(7 


0'  ,  *TAU  =  •  ,F9,6,' 
500)  TAU,C,D 


♦  C  =   '  ,F9Q6,  • 


,AND  D  = 


0'  ,  10F12o6) 
501)  NX 


1*  , 

50  3  ) 

'0','PT  SOURCE  AT  DEPTH •  , F 9„ 6/ ' 0 • ,  •  DENSITIES 

504) XPRI 

0'  ,  «R  =,,F8.>4,hX,1P10E11o3/(»0,,15X,1P10E11o3) 


505)  RMESH(K)  ,  (PT( J,K)  , J  =  l ,NX  ) 

507) 


30(  » 
503) 


•  I,  «R  =«,  F8o4, 


FOLLOWS1 ,30( *  *• ) ) 


505)  RMESH(K) 


(B(J,K)  ,  J=l,NX) 


'0'  ) 

508)  RMFSH  (K) 
509  FORMAT( 1P8E1083) 

WRITE  (7,50<;>  (3(J,K),  J  =  1,NX) 
106  CCNT  INUE 
STOP 
END 
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